package cct.math.polyg;

import java.util.ArrayList;
import java.util.List;

import cct.interfaces.ImplicitFunctionInterface;

/**
 * <p>Title: Molecular Structure Viewer/Editor</p>
 *
 * <p>Description: Computational Chemistry Toolkit</p>
 *
 * <p>Copyright: Copyright (c) 2006</p>
 *
 * <p>Company: ANU</p>
 *
 * @author not attributable
 * @version 1.0
 */
public class Polygonizer {

  // For any edge, if one vertex is inside of the surface and the other is outside of the surface
  //  then the edge intersects the surface
  // For each of the 8 vertices of the cube can be two possible states : either inside or outside of the surface
  // For any cube the are 2^8=256 possible sets of vertex states
  // This table lists the edges intersected by the surface for all 256 possible vertex states
  // There are 12 edges.  For each entry in the table, if edge #n is intersected, then bit #n is set to 1

  static final int aiCubeEdgeFlags[] = { // 256 values
      0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905,
      0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
      0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895,
      0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
      0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35,
      0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
      0x3a0, 0x2a9, 0x1a3, 0x0aa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5,
      0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
      0x460, 0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65,
      0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
      0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5,
      0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
      0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c, 0xe5c, 0xf55,
      0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
      0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5,
      0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
      0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5,
      0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
      0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055,
      0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
      0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5,
      0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
      0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265,
      0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460,
      0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5,
      0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
      0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435,
      0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230,
      0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795,
      0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190,
      0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605,
      0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000
  };

  //  For each of the possible vertex states listed in aiCubeEdgeFlags there is a specific triangulation
//  of the edge intersection points.  a2iTriangleConnectionTable lists all of them in the form of
//  0-5 edge triples with the list terminated by the invalid value -1.
//  For example: a2iTriangleConnectionTable[3] list the 2 triangles formed when corner[0]
//  and corner[1] are inside of the surface, but the rest of the cube is not.
//
//  I found this table in an example program someone wrote long ago.  It was probably generated by hand

  static final int a2iTriangleConnectionTable[][] = { // 256x16
      {
      -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, {
      3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, {
      3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, {
      3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, {
      9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, {
      9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, {
      2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, {
      8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, {
      9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, {
      4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, {
      3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, {
      1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, {
      4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, {
      4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, {
      9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, {
      5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, {
      2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, {
      9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, {
      0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, {
      2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, {
      10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, {
      4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, {
      5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, {
      5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, {
      9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, {
      0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, {
      1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, {
      10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, {
      8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, {
      2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, {
      7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, {
      9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, {
      2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, {
      11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, {
      9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, {
      5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, {
      11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, {
      11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, {
      1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, {
      9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, {
      5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, {
      2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, {
      0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, {
      5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, {
      6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, {
      3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, {
      6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, {
      5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, {
      1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, {
      10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, {
      6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, {
      8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, {
      7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, {
      3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, {
      5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, {
      0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, {
      9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, {
      8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, {
      5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, {
      0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, {
      6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, {
      10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, {
      10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, {
      8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, {
      1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, {
      3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, {
      0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, {
      10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, {
      3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, {
      6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, {
      9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, {
      8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, {
      3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, {
      6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, {
      0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, {
      10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, {
      10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, {
      2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, {
      7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, {
      7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, {
      2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, {
      1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, {
      11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, {
      8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, {
      0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, {
      7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, {
      10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, {
      2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, {
      6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, {
      7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, {
      2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, {
      1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, {
      10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, {
      10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, {
      0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, {
      7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, {
      6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, {
      8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, {
      9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, {
      6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, {
      4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, {
      10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, {
      8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, {
      0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, {
      1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, {
      8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, {
      10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, {
      4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, {
      10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, {
      5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, {
      11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, {
      9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, {
      6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, {
      7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, {
      3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, {
      7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, {
      9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, {
      3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, {
      6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, {
      9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, {
      1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, {
      4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, {
      7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, {
      6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, {
      3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, {
      0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, {
      6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, {
      0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, {
      11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, {
      6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, {
      5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, {
      9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, {
      1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, {
      1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, {
      10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, {
      0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, {
      5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, {
      10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, {
      11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, {
      9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, {
      7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, {
      2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, {
      8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, {
      9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, {
      9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, {
      1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, {
      9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, {
      9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, {
      5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, {
      0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, {
      10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, {
      2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, {
      0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, {
      0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, {
      9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, {
      5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, {
      3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, {
      5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, {
      8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, {
      9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, {
      0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, {
      1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, {
      3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, {
      4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, {
      9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, {
      11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, {
      11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, {
      2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, {
      9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1}, {
      3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1}, {
      1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}, {
      4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1}, {
      4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1}, {
      0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}, {
      3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1}, {
      3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1}, {
      0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}, {
      9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1}, {
      1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {
      -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
  };

  //a2fVertexOffset lists the positions, relative to vertex0, of each of the 8 vertices of a cube
  static final float a2fVertexOffset[][] = { // 8x3
      {
      0.0f, 0.0f, 0.0f}, {
      1.0f, 0.0f, 0.0f}, {
      1.0f, 1.0f, 0.0f}, {
      0.0f, 1.0f, 0.0f}, {
      0.0f, 0.0f, 1.0f}, {
      1.0f, 0.0f, 1.0f}, {
      1.0f, 1.0f, 1.0f}, {
      0.0f, 1.0f, 1.0f}
  };

  //a2iEdgeConnection lists the index of the endpoint vertices for each of the 12 edges of the cube
  static final int a2iEdgeConnection[][] = { // 12x12
      {
      0, 1}, {
      1, 2}, {
      2, 3}, {
      3, 0}, {
      4, 5}, {
      5, 6}, {
      6, 7}, {
      7, 4}, {
      0, 4}, {
      1, 5}, {
      2, 6}, {
      3, 7}
  };

  //a2fEdgeDirection lists the direction vector (vertex1-vertex0) for each edge in the cube
  static final float a2fEdgeDirection[][] = { // 12x3
      {
      1.0f, 0.0f, 0.0f}, {
      0.0f, 1.0f, 0.0f}, {
      -1.0f, 0.0f, 0.0f}, {
      0.0f, -1.0f, 0.0f}, {
      1.0f, 0.0f, 0.0f}, {
      0.0f, 1.0f, 0.0f}, {
      -1.0f, 0.0f, 0.0f}, {
      0.0f, -1.0f, 0.0f}, {
      0.0f, 0.0f, 1.0f}, {
      0.0f, 0.0f, 1.0f}, {
      0.0f, 0.0f, 1.0f}, {
      0.0f, 0.0f, 1.0f}
  };

  //a2iTetrahedronEdgeConnection lists the index of the endpoint vertices for each of the 6 edges of the tetrahedron
  static final int a2iTetrahedronEdgeConnection[][] = { // 6x2
      {
      0, 1}, {
      1, 2}, {
      2, 0}, {
      0, 3}, {
      1, 3}, {
      2, 3}
  };

  //a2iTetrahedronEdgeConnection lists the index of verticies from a cube
  // that made up each of the six tetrahedrons within the cube
  static final int a2iTetrahedronsInACube[][] = { // 6x4
      {
      0, 5, 1, 6}, {
      0, 1, 2, 6}, {
      0, 2, 3, 6}, {
      0, 3, 7, 6}, {
      0, 7, 4, 6}, {
      0, 4, 5, 6},
  };

  // For any edge, if one vertex is inside of the surface and the other is outside of the surface
  //  then the edge intersects the surface
  // For each of the 4 vertices of the tetrahedron can be two possible states : either inside or outside of the surface
  // For any tetrahedron the are 2^4=16 possible sets of vertex states
  // This table lists the edges intersected by the surface for all 16 possible vertex states
  // There are 6 edges.  For each entry in the table, if edge #n is intersected, then bit #n is set to 1
  int aiTetrahedronEdgeFlags[] = { // 16
      0x00, 0x0d, 0x13, 0x1e, 0x26, 0x2b, 0x35, 0x38, 0x38, 0x35, 0x2b, 0x26,
      0x1e, 0x13, 0x0d, 0x00,
  };

  // For each of the possible vertex states listed in aiTetrahedronEdgeFlags there is a specific triangulation
  // of the edge intersection points.  a2iTetrahedronTriangles lists all of them in the form of
  // 0-2 edge triples with the list terminated by the invalid value -1.

  int a2iTetrahedronTriangles[][] = { // 16x7
      {
      -1, -1, -1, -1, -1, -1, -1}, {
      0, 3, 2, -1, -1, -1, -1}, {
      0, 1, 4, -1, -1, -1, -1}, {
      1, 4, 2, 2, 4, 3, -1},

      {
      1, 2, 5, -1, -1, -1, -1}, {
      0, 3, 5, 0, 5, 1, -1}, {
      0, 2, 5, 0, 5, 4, -1}, {
      5, 4, 3, -1, -1, -1, -1},

      {
      3, 4, 5, -1, -1, -1, -1}, {
      4, 5, 0, 5, 2, 0, -1}, {
      1, 5, 0, 5, 3, 0, -1}, {
      5, 2, 1, -1, -1, -1, -1},

      {
      3, 4, 2, 2, 4, 1, -1}, {
      4, 1, 0, -1, -1, -1, -1}, {
      2, 3, 0, -1, -1, -1, -1}, {
      -1, -1, -1, -1, -1, -1, -1},
  };

  protected int[] iDataSetSize = {
      20, 20, 20};
  protected float[] dataOrigin = {
      0, 0, 0};

  private float[] xyzMax = new float[3];

  protected float fStepSize = 0.1f;
  protected float fTargetValue = 1;
  protected ImplicitFunctionInterface function = null;

  protected List Normals = new ArrayList();
  protected List Vertices = new ArrayList();
  protected List Triangles = new ArrayList();

  public Polygonizer() {
    function = new defaultFunction();
  }

  public Polygonizer(float origin[], int dataSetSize[], float stepSize,
                     float isovalue, ImplicitFunctionInterface fun) {
    dataOrigin = origin;
    iDataSetSize = dataSetSize;
    fStepSize = stepSize;
    fTargetValue = isovalue;
    function = fun;

    for (int i = 0; i < 3; i++) {
      xyzMax[i] = dataOrigin[i] + iDataSetSize[i] * stepSize;
    }
  }

  public int noTriangles() {
    if (Triangles == null) {
      return 0;
    }
    return Triangles.size();
  }

  public List getVertices() {
    return Vertices;
  }

  public List getNormals() {
    return Normals;
  }

  public static void main(String[] args) {
    Polygonizer polygonizer = new Polygonizer();
  }

  //vMarchingCubes iterates over the entire dataset, calling vMarchCube on each cube
  public void marchingCubes() {
    int iX, iY, iZ;
    for (iX = 0; iX < iDataSetSize[0]; iX++) {
      for (iY = 0; iY < iDataSetSize[1]; iY++) {
        for (iZ = 0; iZ < iDataSetSize[2]; iZ++) {
          marchCube(dataOrigin[0] + iX * fStepSize,
                    dataOrigin[1] + iY * fStepSize,
                    dataOrigin[2] + iZ * fStepSize, fStepSize);
        }
      }
    }
  }

  public void marchingTetrahedrons() {
    int iX, iY, iZ;
    for (iX = 0; iX < iDataSetSize[0]; iX++) {
      for (iY = 0; iY < iDataSetSize[1]; iY++) {
        for (iZ = 0; iZ < iDataSetSize[2]; iZ++) {
          marchTetrahedrons(dataOrigin[0] + iX * fStepSize,
                            dataOrigin[1] + iY * fStepSize,
                            dataOrigin[2] + iZ * fStepSize, fStepSize);
        }
      }
    }
  }

  //vMarchCube1 performs the Marching Cubes algorithm on a single cube
  void marchCube(float fX, float fY, float fZ, float fScale) {

    int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex,
        iEdgeFlags;
    float fOffset;
    GLvector sColor = new GLvector();
    float[] afCubeValue = new float[8];
    GLvector[] asEdgeVertex = new GLvector[12];
    GLvector[] asEdgeNorm = new GLvector[12];

    for (int i = 0; i < 12; i++) {
      asEdgeVertex[i] = new GLvector();
      asEdgeNorm[i] = new GLvector();
    }

    //Make a local copy of the values at the cube's corners
    for (iVertex = 0; iVertex < 8; iVertex++) {
      afCubeValue[iVertex] = (float) function.eval(fX +
          a2fVertexOffset[iVertex][0] * fScale,
          fY +
          a2fVertexOffset[iVertex][1] * fScale,
          fZ +
          a2fVertexOffset[iVertex][2] * fScale);
    }

    //Find which vertices are inside of the surface and which are outside
    iFlagIndex = 0;
    for (iVertexTest = 0; iVertexTest < 8; iVertexTest++) {
      if (afCubeValue[iVertexTest] <= fTargetValue) {
        iFlagIndex |= 1 << iVertexTest;
      }
    }

    //Find which edges are intersected by the surface
    iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];

    //If the cube is entirely inside or outside of the surface, then there will be no intersections
    if (iEdgeFlags == 0) {
      return;
    }

    //Find the point of intersection of the surface with each edge
    //Then find the normal to the surface at those points
    for (iEdge = 0; iEdge < 12; iEdge++) {
      //if there is an intersection on this edge
      if ( (iEdgeFlags & (1 << iEdge)) != 0) {
        fOffset = fGetOffset(afCubeValue[a2iEdgeConnection[iEdge][0]],
                             afCubeValue[a2iEdgeConnection[iEdge][1]],
                             fTargetValue);

        asEdgeVertex[iEdge].fX = fX +
            (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][0] +
             fOffset * a2fEdgeDirection[iEdge][0]) * fScale;
        asEdgeVertex[iEdge].fY = fY +
            (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][1] +
             fOffset * a2fEdgeDirection[iEdge][1]) * fScale;
        asEdgeVertex[iEdge].fZ = fZ +
            (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][2] +
             fOffset * a2fEdgeDirection[iEdge][2]) * fScale;

        getNormal(asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX,
                  asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
      }
    }

    //Draw the triangles that were found.  There can be up to five per cube
    for (iTriangle = 0; iTriangle < 5; iTriangle++) {
      if (a2iTriangleConnectionTable[iFlagIndex][3 * iTriangle] < 0) {
        break;
      }

      Triangle triangle = new Triangle();

      for (iCorner = 0; iCorner < 3; iCorner++) {
        iVertex = a2iTriangleConnectionTable[iFlagIndex][3 * iTriangle +
            iCorner];

        //vGetColor(sColor, asEdgeVertex[iVertex], asEdgeNorm[iVertex]);
        //glColor3f(sColor.fX, sColor.fY, sColor.fZ);

        //glNormal3f(asEdgeNorm[iVertex].fX, asEdgeNorm[iVertex].fY,
        //           asEdgeNorm[iVertex].fZ);
        GLvector normal = new GLvector(asEdgeNorm[iVertex].fX,
                                       asEdgeNorm[iVertex].fY,
                                       asEdgeNorm[iVertex].fZ);
        Normals.add(normal);

        //glVertex3f(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY,
        //           asEdgeVertex[iVertex].fZ);

        GLvector vertex = new GLvector(asEdgeVertex[iVertex].fX,
                                       asEdgeVertex[iVertex].fY,
                                       asEdgeVertex[iVertex].fZ);
        Vertices.add(vertex);

        triangle.addVertex(vertex, normal);
      }

      Triangles.add(triangle);
    }
  }

  //vMarchingCubes iterates over the entire dataset, calling vMarchCube on each cube
  public void marchingCubes2() {
    int iX, iY, iZ;
    Cube[][] slab_0 = new Cube[iDataSetSize[0]][iDataSetSize[1]];
    Cube[][] slab_1 = new Cube[iDataSetSize[0]][iDataSetSize[1]];

    for (iX = 0; iX < iDataSetSize[0]; iX++) {
      for (iY = 0; iY < iDataSetSize[1]; iY++) {
        slab_0[iX][iY] = new Cube();
        slab_0[iX][iY].setStep(fStepSize);
        slab_1[iX][iY] = new Cube();
        slab_1[iX][iY].setStep(fStepSize);
      }
    }

    for (iZ = 0; iZ < iDataSetSize[2]; iZ++) {
      float fZ = dataOrigin[2] + iZ * fStepSize;

      for (iX = 0; iX < iDataSetSize[0]; iX++) {
        float fX = dataOrigin[0] + iX * fStepSize;
        for (iY = 0; iY < iDataSetSize[1]; iY++) {
          float fY = dataOrigin[1] + iY * fStepSize;

          Cube cube = slab_1[iX][iY];
          cube.setOrigin(fX, fY, fZ);

          setupCube(slab_0, slab_1, iX, iY, iZ);

        }
      }
      // --- Switch slabs
      Cube[][] temp = slab_0;
      slab_0 = slab_1;
      slab_1 = temp;
    }
  }

  void setupCube(Cube[][] slab_0, Cube[][] slab_1, int iX, int iY, int iZ) {

    // --- Set function values

    Cube cube = slab_1[iX][iY];

    for (int iVertex = 0; iVertex < 8; iVertex++) {
      cube.afCubeValue[iVertex] = (float) function.eval(cube.originX +
          a2fVertexOffset[iVertex][0] * cube.step,
          cube.originY +
          a2fVertexOffset[iVertex][1] * cube.step,
          cube.originZ +
          a2fVertexOffset[iVertex][2] * cube.step);
    }

    if (true) {

      if (iX == 0 && iY == 0) {
        // --- Set 4,5, 6, 7th vertices
        cube.afCubeValue[4] = (float) function.eval(cube.originX + a2fVertexOffset[4][0] * cube.step,
            cube.originY + a2fVertexOffset[4][1] * cube.step,
            cube.originZ + a2fVertexOffset[4][2] * cube.step);
        cube.afCubeValue[5] = (float) function.eval(cube.originX + a2fVertexOffset[5][0] * cube.step,
            cube.originY + a2fVertexOffset[5][1] * cube.step,
            cube.originZ + a2fVertexOffset[5][2] * cube.step);
        cube.afCubeValue[6] = (float) function.eval(cube.originX + a2fVertexOffset[6][0] * cube.step,
            cube.originY + a2fVertexOffset[6][1] * cube.step,
            cube.originZ + a2fVertexOffset[6][2] * cube.step);
        cube.afCubeValue[7] = (float) function.eval(cube.originX + a2fVertexOffset[7][0] * cube.step,
            cube.originY + a2fVertexOffset[7][1] * cube.step,
            cube.originZ + a2fVertexOffset[7][2] * cube.step);
      }
      else if (iX == 0) {

        Cube cube_0 = slab_1[iX][iY - 1];
        cube.afCubeValue[4] = cube_0.afCubeValue[7];
        cube.afCubeValue[5] = cube_0.afCubeValue[6];
        cube.afCubeValue[6] = (float) function.eval(cube.originX + a2fVertexOffset[6][0] * fStepSize,
            cube.originY + a2fVertexOffset[6][1] * fStepSize,
            cube.originZ + a2fVertexOffset[6][2] * fStepSize);
        cube.afCubeValue[7] = (float) function.eval(cube.originX + a2fVertexOffset[7][0] * fStepSize,
            cube.originY + a2fVertexOffset[7][1] * fStepSize,
            cube.originZ + a2fVertexOffset[7][2] * fStepSize);
      }
      else if (iY == 0) {

        Cube cube_0 = slab_1[iX - 1][iY];
        cube.afCubeValue[4] = cube_0.afCubeValue[5];
        cube.afCubeValue[7] = cube_0.afCubeValue[6];
        cube.afCubeValue[5] = (float) function.eval(cube.originX + a2fVertexOffset[5][0] * fStepSize,
            cube.originY + a2fVertexOffset[5][1] * fStepSize,
            cube.originZ + a2fVertexOffset[5][2] * fStepSize);
        cube.afCubeValue[6] = (float) function.eval(cube.originX + a2fVertexOffset[6][0] * fStepSize,
            cube.originY + a2fVertexOffset[6][1] * fStepSize,
            cube.originZ + a2fVertexOffset[6][2] * fStepSize);
      }
      else {

        Cube cube_0 = slab_1[iX][iY - 1];
        cube.afCubeValue[4] = cube_0.afCubeValue[7];
        cube.afCubeValue[5] = cube_0.afCubeValue[6];
        cube_0 = slab_1[iX - 1][iY];
        cube.afCubeValue[7] = cube_0.afCubeValue[6];
        cube.afCubeValue[6] = (float) function.eval(cube.originX + a2fVertexOffset[6][0] * fStepSize,
            cube.originY + a2fVertexOffset[6][1] * fStepSize,
            cube.originZ + a2fVertexOffset[6][2] * fStepSize);
      }

      // --- Process slab
      if (iZ == 0) {
        if (iX == 0 && iY == 0) {
          // --- Set 4,5, 6, 7th vertices
          cube.afCubeValue[0] = (float) function.eval(cube.originX + a2fVertexOffset[0][0] * fStepSize,
              cube.originY + a2fVertexOffset[0][1] * fStepSize,
              cube.originZ + a2fVertexOffset[0][2] * fStepSize);
          cube.afCubeValue[1] = (float) function.eval(cube.originX + a2fVertexOffset[1][0] * fStepSize,
              cube.originY + a2fVertexOffset[1][1] * fStepSize,
              cube.originZ + a2fVertexOffset[1][2] * fStepSize);
          cube.afCubeValue[2] = (float) function.eval(cube.originX + a2fVertexOffset[2][0] * fStepSize,
              cube.originY + a2fVertexOffset[2][1] * fStepSize,
              cube.originZ + a2fVertexOffset[2][2] * fStepSize);
          cube.afCubeValue[3] = (float) function.eval(cube.originX + a2fVertexOffset[3][0] * fStepSize,
              cube.originY + a2fVertexOffset[3][1] * fStepSize,
              cube.originZ + a2fVertexOffset[3][2] * fStepSize);
        }
        else if (iX == 0) {
          Cube cube_0 = slab_1[iX][iY - 1];
          cube.afCubeValue[0] = cube_0.afCubeValue[3];
          cube.afCubeValue[1] = cube_0.afCubeValue[2];
          cube.afCubeValue[2] = (float) function.eval(cube.originX + a2fVertexOffset[2][0] * fStepSize,
              cube.originY + a2fVertexOffset[2][1] * fStepSize,
              cube.originZ + a2fVertexOffset[2][2] * fStepSize);
          cube.afCubeValue[3] = (float) function.eval(cube.originX + a2fVertexOffset[3][0] * fStepSize,
              cube.originY + a2fVertexOffset[3][1] * fStepSize,
              cube.originZ + a2fVertexOffset[3][2] * fStepSize);
        }
        else if (iY == 0) {
          Cube cube_0 = slab_1[iX - 1][iY];
          cube.afCubeValue[0] = cube_0.afCubeValue[1];
          cube.afCubeValue[3] = cube_0.afCubeValue[2];
          cube.afCubeValue[1] = (float) function.eval(cube.originX + a2fVertexOffset[1][0] * fStepSize,
              cube.originY + a2fVertexOffset[1][1] * fStepSize,
              cube.originZ + a2fVertexOffset[1][2] * fStepSize);
          cube.afCubeValue[2] = (float) function.eval(cube.originX + a2fVertexOffset[2][0] * fStepSize,
              cube.originY + a2fVertexOffset[2][1] * fStepSize,
              cube.originZ + a2fVertexOffset[2][2] * fStepSize);
        }
        else {
          Cube cube_0 = slab_1[iX][iY - 1];
          cube.afCubeValue[0] = cube_0.afCubeValue[3];
          cube.afCubeValue[1] = cube_0.afCubeValue[2];
          cube_0 = slab_1[iX - 1][iY];
          cube.afCubeValue[3] = cube_0.afCubeValue[2];
          cube.afCubeValue[2] = (float) function.eval(cube.originX + a2fVertexOffset[2][0] * fStepSize,
              cube.originY + a2fVertexOffset[2][1] * fStepSize,
              cube.originZ + a2fVertexOffset[2][2] * fStepSize);
        }

      }

      else {
        Cube cube_0 = slab_0[iX][iY];
        cube.afCubeValue[0] = cube_0.afCubeValue[4];
        cube.afCubeValue[1] = cube_0.afCubeValue[5];
        cube.afCubeValue[2] = cube_0.afCubeValue[6];
        cube.afCubeValue[3] = cube_0.afCubeValue[7];
      }
    }

    //Find which vertices are inside of the surface and which are outside
    int iEdgeFlags = cube.getEdgeFlags(fTargetValue);

    //If the cube is entirely inside or outside of the surface, then there will be no intersections
    if (iEdgeFlags == 0) {
      return;
    }

    //Find the point of intersection of the surface with each edge
    //Then find the normal to the surface at those points

    for (int iEdge = 0; iEdge < 12; iEdge++) {

      //if there is an intersection on this edge
      if ( (iEdgeFlags & (1 << iEdge)) != 0) {
        float fOffset = fGetOffset(cube.afCubeValue[a2iEdgeConnection[iEdge][0]],
                                   cube.afCubeValue[a2iEdgeConnection[iEdge][1]],
                                   fTargetValue);

        cube.asEdgeVertex[iEdge].fX = cube.originX + (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][0] +
            fOffset * a2fEdgeDirection[iEdge][0]) * fStepSize;
        cube.asEdgeVertex[iEdge].fY = cube.originY + (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][1] +
            fOffset * a2fEdgeDirection[iEdge][1]) * fStepSize;
        cube.asEdgeVertex[iEdge].fZ = cube.originZ + (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][2] +
            fOffset * a2fEdgeDirection[iEdge][2]) * fStepSize;

        getNormal(cube.asEdgeNorm[iEdge], cube.asEdgeVertex[iEdge].fX,
                  cube.asEdgeVertex[iEdge].fY, cube.asEdgeVertex[iEdge].fZ);
      }
    }

    if (false) {
      // More smart stuff, if any

      if (iX == 0 && iY == 0) {
        if ( (iEdgeFlags & (1 << 4)) != 0) {
          setupEdge(cube, 4);
        }
        if ( (iEdgeFlags & (1 << 5)) != 0) {
          setupEdge(cube, 5);
        }
        if ( (iEdgeFlags & (1 << 6)) != 0) {
          setupEdge(cube, 6);
        }
        if ( (iEdgeFlags & (1 << 7)) != 0) {
          setupEdge(cube, 7);
        }
        if ( (iEdgeFlags & (1 << 8)) != 0) {
          setupEdge(cube, 8);
        }
        if ( (iEdgeFlags & (1 << 9)) != 0) {
          setupEdge(cube, 9);
        }
        if ( (iEdgeFlags & (1 << 10)) != 0) {
          setupEdge(cube, 10);
        }
        if ( (iEdgeFlags & (1 << 11)) != 0) {
          setupEdge(cube, 11);
        }
      }

      else if (iX == 0) {

        if ( (iEdgeFlags & (1 << 5)) != 0) {
          setupEdge(cube, 5);
        }
        if ( (iEdgeFlags & (1 << 6)) != 0) {
          setupEdge(cube, 6);
        }
        if ( (iEdgeFlags & (1 << 7)) != 0) {
          setupEdge(cube, 7);
        }
        if ( (iEdgeFlags & (1 << 10)) != 0) {
          setupEdge(cube, 10);
        }
        if ( (iEdgeFlags & (1 << 11)) != 0) {
          setupEdge(cube, 11);
        }

        Cube cube_0 = slab_1[iX][iY - 1];

        if ( (iEdgeFlags & (1 << 4)) != 0) {
          setupEdge(cube, 4, cube_0, 6);
        }
        if ( (iEdgeFlags & (1 << 8)) != 0) {
          setupEdge(cube, 8, cube_0, 11);
        }
        if ( (iEdgeFlags & (1 << 9)) != 0) {
          setupEdge(cube, 9, cube_0, 10);
        }
      }

      else if (iY == 0) {

        if ( (iEdgeFlags & (1 << 4)) != 0) {
          setupEdge(cube, 4);
        }
        if ( (iEdgeFlags & (1 << 5)) != 0) {
          setupEdge(cube, 5);
        }
        if ( (iEdgeFlags & (1 << 6)) != 0) {
          setupEdge(cube, 6);
        }
        if ( (iEdgeFlags & (1 << 9)) != 0) {
          setupEdge(cube, 9);
        }
        if ( (iEdgeFlags & (1 << 10)) != 0) {
          setupEdge(cube, 10);
        }

        Cube cube_0 = slab_1[iX - 1][iY];

        if ( (iEdgeFlags & (1 << 7)) != 0) {
          setupEdge(cube, 7, cube_0, 5);
        }
        if ( (iEdgeFlags & (1 << 8)) != 0) {
          setupEdge(cube, 8, cube_0, 9);
        }
        if ( (iEdgeFlags & (1 << 11)) != 0) {
          setupEdge(cube, 11, cube_0, 10);
        }
      }

      else {
        if ( (iEdgeFlags & (1 << 5)) != 0) {
          setupEdge(cube, 5);
        }
        if ( (iEdgeFlags & (1 << 6)) != 0) {
          setupEdge(cube, 6);
        }
        if ( (iEdgeFlags & (1 << 10)) != 0) {
          setupEdge(cube, 10);
        }

        Cube cube_0 = slab_1[iX][iY - 1];

        if ( (iEdgeFlags & (1 << 4)) != 0) {
          setupEdge(cube, 4, cube_0, 6);
        }
        if ( (iEdgeFlags & (1 << 8)) != 0) {
          setupEdge(cube, 8, cube_0, 11);
        }
        if ( (iEdgeFlags & (1 << 9)) != 0) {
          setupEdge(cube, 9, cube_0, 10);
        }

        cube_0 = slab_1[iX - 1][iY];

        if ( (iEdgeFlags & (1 << 7)) != 0) {
          setupEdge(cube, 7, cube_0, 5);
        }
        if ( (iEdgeFlags & (1 << 11)) != 0) {
          setupEdge(cube, 11, cube_0, 10);
        }
      }

      // --- Process slab
      if (iZ == 0) {
        if (iX == 0 && iY == 0) {

          if ( (iEdgeFlags & (1 << 0)) != 0) {
            setupEdge(cube, 0);
          }
          if ( (iEdgeFlags & (1 << 1)) != 0) {
            setupEdge(cube, 1);
          }
          if ( (iEdgeFlags & (1 << 2)) != 0) {
            setupEdge(cube, 2);
          }
          if ( (iEdgeFlags & (1 << 3)) != 0) {
            setupEdge(cube, 3);
          }
        }

        else if (iX == 0) {

          if ( (iEdgeFlags & (1 << 1)) != 0) {
            setupEdge(cube, 1);
          }
          if ( (iEdgeFlags & (1 << 2)) != 0) {
            setupEdge(cube, 2);
          }
          if ( (iEdgeFlags & (1 << 3)) != 0) {
            setupEdge(cube, 3);
          }

          Cube cube_0 = slab_0[iX][iY - 1];

          if ( (iEdgeFlags & (1 << 0)) != 0) {
            setupEdge(cube, 0, cube_0, 2);
          }

        }
        else if (iY == 0) {

          if ( (iEdgeFlags & (1 << 0)) != 0) {
            setupEdge(cube, 0);
          }
          if ( (iEdgeFlags & (1 << 1)) != 0) {
            setupEdge(cube, 1);
          }
          if ( (iEdgeFlags & (1 << 2)) != 0) {
            setupEdge(cube, 2);
          }

          Cube cube_0 = slab_0[iX - 1][iY];

          if ( (iEdgeFlags & (1 << 3)) != 0) {
            setupEdge(cube, 3, cube_0, 1);
          }
        }

        else {
          if ( (iEdgeFlags & (1 << 1)) != 0) {
            setupEdge(cube, 1);
          }
          if ( (iEdgeFlags & (1 << 2)) != 0) {
            setupEdge(cube, 2);
          }

          Cube cube_0 = slab_0[iX][iY - 1];

          if ( (iEdgeFlags & (1 << 0)) != 0) {
            setupEdge(cube, 0, cube_0, 2);
          }

          cube_0 = slab_0[iX - 1][iY];

          if ( (iEdgeFlags & (1 << 3)) != 0) {
            setupEdge(cube, 3, cube_0, 1);
          }
        }
      }

      else {
        Cube cube_0 = slab_0[iX][iY];
        if ( (iEdgeFlags & (1 << 0)) != 0) {
          setupEdge(cube, 0, cube_0, 4);
        }
        if ( (iEdgeFlags & (1 << 1)) != 0) {
          setupEdge(cube, 1, cube_0, 5);
        }
        if ( (iEdgeFlags & (1 << 2)) != 0) {
          setupEdge(cube, 2, cube_0, 6);
        }
        if ( (iEdgeFlags & (1 << 3)) != 0) {
          setupEdge(cube, 3, cube_0, 7);
        }

      }
    }
    //Draw the triangles that were found.  There can be up to five per cube
    int iFlagIndex = cube.getFlagIndex(fTargetValue);
    for (int iTriangle = 0; iTriangle < 5; iTriangle++) {
      if (a2iTriangleConnectionTable[iFlagIndex][3 * iTriangle] < 0) {
        break;
      }

      Triangle triangle = new Triangle();

      for (int iCorner = 0; iCorner < 3; iCorner++) {
        int iVertex = a2iTriangleConnectionTable[iFlagIndex][3 * iTriangle + iCorner];

        GLvector normal = new GLvector(cube.asEdgeNorm[iVertex].fX,
                                       cube.asEdgeNorm[iVertex].fY,
                                       cube.asEdgeNorm[iVertex].fZ);
        Normals.add(normal);

        GLvector vertex = new GLvector(cube.asEdgeVertex[iVertex].fX,
                                       cube.asEdgeVertex[iVertex].fY,
                                       cube.asEdgeVertex[iVertex].fZ);
        Vertices.add(vertex);

        triangle.addVertex(vertex, normal);
      }

      Triangles.add(triangle);
    }

  }

  void setupEdge(Cube cube, int iEdge) {
    float fOffset = fGetOffset(cube.afCubeValue[a2iEdgeConnection[iEdge][0]], cube.afCubeValue[a2iEdgeConnection[iEdge][1]],
                               fTargetValue);

    cube.asEdgeVertex[iEdge].fX = cube.originX + (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][0] +
                                                  fOffset * a2fEdgeDirection[iEdge][0]) * fStepSize;
    cube.asEdgeVertex[iEdge].fY = cube.originY + (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][1] +
                                                  fOffset * a2fEdgeDirection[iEdge][1]) * fStepSize;
    cube.asEdgeVertex[iEdge].fZ = cube.originZ + (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][2] +
                                                  fOffset * a2fEdgeDirection[iEdge][2]) * fStepSize;

    getNormal(cube.asEdgeNorm[iEdge], cube.asEdgeVertex[iEdge].fX, cube.asEdgeVertex[iEdge].fY, cube.asEdgeVertex[iEdge].fZ);
  }

  void setupEdge(Cube cube, int iEdge, Cube cube_2, int iEdge_2) {
    cube.asEdgeVertex[iEdge].fX = cube_2.asEdgeVertex[iEdge_2].fX;
    cube.asEdgeVertex[iEdge].fY = cube_2.asEdgeVertex[iEdge_2].fY;
    cube.asEdgeVertex[iEdge].fZ = cube_2.asEdgeVertex[iEdge_2].fZ;

    cube.asEdgeNorm[iEdge].fX = cube_2.asEdgeNorm[iEdge_2].fX;
    cube.asEdgeNorm[iEdge].fY = cube_2.asEdgeNorm[iEdge_2].fY;
    cube.asEdgeNorm[iEdge].fZ = cube_2.asEdgeNorm[iEdge_2].fZ;
  }

  void marchCube2(Cube cube, float fX, float fY, float fZ, float fScale) {

    int iCorner, iVertex, iVertexTest, iEdge, iTriangle, iFlagIndex, iEdgeFlags;
    float fOffset;
    GLvector sColor = new GLvector();
    float[] afCubeValue = new float[8];
    GLvector[] asEdgeVertex = new GLvector[12];
    GLvector[] asEdgeNorm = new GLvector[12];

    for (int i = 0; i < 12; i++) {
      asEdgeVertex[i] = new GLvector();
      asEdgeNorm[i] = new GLvector();
    }

    //Make a local copy of the values at the cube's corners
    /*
     0        0.0f, 0.0f, 0.0f}, {
     1      1.0f, 0.0f, 0.0f}, {
     2      1.0f, 1.0f, 0.0f}, {
     3      0.0f, 1.0f, 0.0f}, {
     4      0.0f, 0.0f, 1.0f}, {
     5      1.0f, 0.0f, 1.0f}, {
     6      1.0f, 1.0f, 1.0f}, {
     7      0.0f, 1.0f, 1.0f}
     */
    for (iVertex = 0; iVertex < 8; iVertex++) {
      afCubeValue[iVertex] = (float) function.eval(fX +
          a2fVertexOffset[iVertex][0] * fScale,
          fY +
          a2fVertexOffset[iVertex][1] * fScale,
          fZ +
          a2fVertexOffset[iVertex][2] * fScale);
    }

    //Find which vertices are inside of the surface and which are outside
    iFlagIndex = 0;
    for (iVertexTest = 0; iVertexTest < 8; iVertexTest++) {
      if (afCubeValue[iVertexTest] <= fTargetValue) {
        iFlagIndex |= 1 << iVertexTest;
      }
    }

    //Find which edges are intersected by the surface
    iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];

    //If the cube is entirely inside or outside of the surface, then there will be no intersections
    if (iEdgeFlags == 0) {
      return;
    }

    //Find the point of intersection of the surface with each edge
    //Then find the normal to the surface at those points
    for (iEdge = 0; iEdge < 12; iEdge++) {
      //if there is an intersection on this edge
      if ( (iEdgeFlags & (1 << iEdge)) != 0) {
        fOffset = fGetOffset(afCubeValue[a2iEdgeConnection[iEdge][0]],
                             afCubeValue[a2iEdgeConnection[iEdge][1]],
                             fTargetValue);

        asEdgeVertex[iEdge].fX = fX +
            (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][0] +
             fOffset * a2fEdgeDirection[iEdge][0]) * fScale;
        asEdgeVertex[iEdge].fY = fY +
            (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][1] +
             fOffset * a2fEdgeDirection[iEdge][1]) * fScale;
        asEdgeVertex[iEdge].fZ = fZ +
            (a2fVertexOffset[a2iEdgeConnection[iEdge][0]][2] +
             fOffset * a2fEdgeDirection[iEdge][2]) * fScale;

        getNormal(asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX,
                  asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
      }
    }

    //Draw the triangles that were found.  There can be up to five per cube
    for (iTriangle = 0; iTriangle < 5; iTriangle++) {
      if (a2iTriangleConnectionTable[iFlagIndex][3 * iTriangle] < 0) {
        break;
      }

      Triangle triangle = new Triangle();

      for (iCorner = 0; iCorner < 3; iCorner++) {
        iVertex = a2iTriangleConnectionTable[iFlagIndex][3 * iTriangle +
            iCorner];

        GLvector normal = new GLvector(asEdgeNorm[iVertex].fX,
                                       asEdgeNorm[iVertex].fY,
                                       asEdgeNorm[iVertex].fZ);
        Normals.add(normal);

        GLvector vertex = new GLvector(asEdgeVertex[iVertex].fX,
                                       asEdgeVertex[iVertex].fY,
                                       asEdgeVertex[iVertex].fZ);
        Vertices.add(vertex);

        triangle.addVertex(vertex, normal);
      }

      Triangles.add(triangle);
    }
  }

//vGetNormal() finds the gradient of the scalar field at a point
//This gradient can be used as a very accurate vertx normal for lighting calculations
  void getNormal(GLvector rfNormal, float fX, float fY, float fZ) {
    float delta = 0.01f;

    if (fX - delta < this.dataOrigin[0]) {
      rfNormal.fX = (float) (function.eval(fX, fY, fZ) - function.eval(fX + delta, fY, fZ)) / delta;
    }
    else if (fX + delta > xyzMax[0]) {
      rfNormal.fX = (float) (function.eval(fX - delta, fY, fZ) - function.eval(fX, fY, fZ)) / delta;
    }
    else {
      rfNormal.fX = (float) (function.eval(fX - delta, fY, fZ) - function.eval(fX + delta, fY, fZ)) / (2 * delta);
    }

    if (fY - delta < this.dataOrigin[1]) {
      rfNormal.fY = (float) (function.eval(fX, fY, fZ) - function.eval(fX, fY + delta, fZ)) / delta;
    }
    else if (fY + delta > xyzMax[1]) {
      rfNormal.fY = (float) (function.eval(fX, fY - delta, fZ) - function.eval(fX, fY, fZ)) / delta;
    }
    else {
      rfNormal.fY = (float) (function.eval(fX, fY - delta, fZ) - function.eval(fX, fY + delta, fZ)) / (2 * delta);
    }

    if (fZ - delta < this.dataOrigin[2]) {
      rfNormal.fZ = (float) (function.eval(fX, fY, fZ) - function.eval(fX, fY, fZ + delta)) / delta;
    }
    else if (fZ + delta > xyzMax[2]) {
      rfNormal.fZ = (float) (function.eval(fX, fY, fZ - delta) - function.eval(fX, fY, fZ)) / delta;
    }
    else {
      rfNormal.fZ = (float) (function.eval(fX, fY, fZ - delta) - function.eval(fX, fY, fZ + delta)) / (2 * delta);
    }

    /*
         rfNormal.fX = function.eval(fX - 0.01f, fY, fZ) -        function.eval(fX + 0.01f, fY, fZ);
         rfNormal.fY = function.eval(fX, fY - 0.01f, fZ) -        function.eval(fX, fY + 0.01f, fZ);
         rfNormal.fZ = function.eval(fX, fY, fZ - 0.01f) -        function.eval(fX, fY, fZ + 0.01f);
     */

    normalizeVector(rfNormal, rfNormal);
  }

  void normalizeVector(GLvector rfVectorResult, GLvector rfVectorSource) {
    float fOldLength;
    float fScale;

    fOldLength = (float) Math.sqrt( (rfVectorSource.fX * rfVectorSource.fX) +
                                   (rfVectorSource.fY * rfVectorSource.fY) +
                                   (rfVectorSource.fZ * rfVectorSource.fZ));

    if (fOldLength == 0.0) {
      rfVectorResult.fX = rfVectorSource.fX;
      rfVectorResult.fY = rfVectorSource.fY;
      rfVectorResult.fZ = rfVectorSource.fZ;
    }
    else {
      fScale = 1.0f / fOldLength;
      rfVectorResult.fX = rfVectorSource.fX * fScale;
      rfVectorResult.fY = rfVectorSource.fY * fScale;
      rfVectorResult.fZ = rfVectorSource.fZ * fScale;
    }
  }

//fGetOffset finds the approximate point of intersection of the surface
// between two points with the values fValue1 and fValue2
  float fGetOffset(float fValue1, float fValue2, float fValueDesired) {
    //double fDelta = fValue2 - fValue1;
    float fDelta = fValue2 - fValue1;

    if (fDelta == 0.0f) {
      return 0.5f;
    }
    return (fValueDesired - fValue1) / fDelta;
  }

//vMarchCube2 performs the Marching Tetrahedrons algorithm on a single cube by making six calls to vMarchTetrahedron
  public void marchTetrahedrons(float fX, float fY, float fZ, float fScale) {
    int iVertex, iTetrahedron, iVertexInACube;
    GLvector[] asCubePosition = new GLvector[8];
    float[] afCubeValue = new float[8];
    GLvector[] asTetrahedronPosition = new GLvector[4];
    float[] afTetrahedronValue = new float[4];

    for (int i = 0; i < 8; i++) {
      asCubePosition[i] = new GLvector();
    }
    for (int i = 0; i < 4; i++) {
      asTetrahedronPosition[i] = new GLvector();
    }

    //Make a local copy of the cube's corner positions
    for (iVertex = 0; iVertex < 8; iVertex++) {
      asCubePosition[iVertex].fX = fX + a2fVertexOffset[iVertex][0] * fScale;
      asCubePosition[iVertex].fY = fY + a2fVertexOffset[iVertex][1] * fScale;
      asCubePosition[iVertex].fZ = fZ + a2fVertexOffset[iVertex][2] * fScale;
    }

    //Make a local copy of the cube's corner values
    for (iVertex = 0; iVertex < 8; iVertex++) {
      afCubeValue[iVertex] = (float) function.eval(asCubePosition[iVertex].fX,
          asCubePosition[iVertex].fY,
          asCubePosition[iVertex].fZ);
    }

    for (iTetrahedron = 0; iTetrahedron < 6; iTetrahedron++) {
      for (iVertex = 0; iVertex < 4; iVertex++) {
        iVertexInACube = a2iTetrahedronsInACube[iTetrahedron][iVertex];
        asTetrahedronPosition[iVertex].fX = asCubePosition[iVertexInACube].fX;
        asTetrahedronPosition[iVertex].fY = asCubePosition[iVertexInACube].fY;
        asTetrahedronPosition[iVertex].fZ = asCubePosition[iVertexInACube].fZ;
        afTetrahedronValue[iVertex] = afCubeValue[iVertexInACube];
      }
      marchTetrahedron(asTetrahedronPosition, afTetrahedronValue);
    }
  }

//marchTetrahedron performs the Marching Tetrahedrons algorithm on a single tetrahedron
  void marchTetrahedron(GLvector pasTetrahedronPosition[],
                        float pafTetrahedronValue[]) {
    //extern GLint aiTetrahedronEdgeFlags[16];
    //extern GLint a2iTetrahedronTriangles[16][7];

    int iEdge, iVert0, iVert1, iEdgeFlags, iTriangle, iCorner, iVertex,
        iFlagIndex = 0;
    float fOffset, fInvOffset, fValue = 0.0f;
    GLvector[] asEdgeVertex = new GLvector[6];
    GLvector[] asEdgeNorm = new GLvector[6];
    GLvector sColor;

    for (int i = 0; i < 6; i++) {
      asEdgeVertex[i] = new GLvector();
      asEdgeNorm[i] = new GLvector();
    }

    //Find which vertices are inside of the surface and which are outside
    for (iVertex = 0; iVertex < 4; iVertex++) {
      if (pafTetrahedronValue[iVertex] <= fTargetValue) {
        iFlagIndex |= 1 << iVertex;
      }
    }

    //Find which edges are intersected by the surface
    iEdgeFlags = aiTetrahedronEdgeFlags[iFlagIndex];

    //If the tetrahedron is entirely inside or outside of the surface, then there will be no intersections
    if (iEdgeFlags == 0) {
      return;
    }
    //Find the point of intersection of the surface with each edge
    // Then find the normal to the surface at those points
    for (iEdge = 0; iEdge < 6; iEdge++) {
      //if there is an intersection on this edge
      if ( (iEdgeFlags & (1 << iEdge)) != 0) {
        iVert0 = a2iTetrahedronEdgeConnection[iEdge][0];
        iVert1 = a2iTetrahedronEdgeConnection[iEdge][1];
        fOffset = fGetOffset(pafTetrahedronValue[iVert0],
                             pafTetrahedronValue[iVert1], fTargetValue);
        fInvOffset = 1.0f - fOffset;

        asEdgeVertex[iEdge].fX = fInvOffset * pasTetrahedronPosition[iVert0].fX +
            fOffset * pasTetrahedronPosition[iVert1].fX;
        asEdgeVertex[iEdge].fY = fInvOffset * pasTetrahedronPosition[iVert0].fY +
            fOffset * pasTetrahedronPosition[iVert1].fY;
        asEdgeVertex[iEdge].fZ = fInvOffset * pasTetrahedronPosition[iVert0].fZ +
            fOffset * pasTetrahedronPosition[iVert1].fZ;

        getNormal(asEdgeNorm[iEdge], asEdgeVertex[iEdge].fX,
                  asEdgeVertex[iEdge].fY, asEdgeVertex[iEdge].fZ);
      }
    }
    //Draw the triangles that were found.  There can be up to 2 per tetrahedron
    for (iTriangle = 0; iTriangle < 2; iTriangle++) {
      if (a2iTetrahedronTriangles[iFlagIndex][3 * iTriangle] < 0) {
        break;
      }

      Triangle triangle = new Triangle();

      for (iCorner = 0; iCorner < 3; iCorner++) {
        iVertex = a2iTetrahedronTriangles[iFlagIndex][3 * iTriangle + iCorner];

        //vGetColor(sColor, asEdgeVertex[iVertex], asEdgeNorm[iVertex]);
        //glColor3f(sColor.fX, sColor.fY, sColor.fZ);

        //glNormal3f(asEdgeNorm[iVertex].fX, asEdgeNorm[iVertex].fY,
        //           asEdgeNorm[iVertex].fZ);
        GLvector normal = new GLvector(asEdgeNorm[iVertex].fX,
                                       asEdgeNorm[iVertex].fY,
                                       asEdgeNorm[iVertex].fZ);
        Normals.add(normal);

        //glVertex3f(asEdgeVertex[iVertex].fX, asEdgeVertex[iVertex].fY,
        //           asEdgeVertex[iVertex].fZ);
        GLvector vertex = new GLvector(asEdgeVertex[iVertex].fX,
                                       asEdgeVertex[iVertex].fY,
                                       asEdgeVertex[iVertex].fZ);
        Vertices.add(vertex);

        triangle.addVertex(vertex, normal);

      }

      Triangles.add(triangle);
    }
  }

  public void dispose() {
    Normals.clear();
    Vertices.clear();
    Triangles.clear();
    Normals = null;
    Vertices = null;
    Triangles = null;
  }

  private class Triangle {
    GLvector[] vertices = new GLvector[3];
    GLvector[] normals = new GLvector[3];
    int index = 0;

    public void addVertex(GLvector vertex, GLvector normal) {
      vertices[index] = vertex;
      normals[index] = normal;
      ++index;
    }
  }

  private class Cube {
    float originX, originY, originZ, step;
    float[] afCubeValue = new float[8];
    GLvector[] asEdgeVertex = new GLvector[12];
    GLvector[] asEdgeNorm = new GLvector[12];

    public Cube() {
      for (int i = 0; i < 12; i++) {
        asEdgeVertex[i] = new GLvector();
        asEdgeNorm[i] = new GLvector();
      }
    }

    public void setStep(float st) {
      step = st;
    }

    public void setOrigin(float x, float y, float z) {
      originX = x;
      originY = y;
      originZ = z;
    }

    public int getFlagIndex(float fTargetValue) {
      //Find which vertices are inside of the surface and which are outside
      int iFlagIndex = 0;
      for (int iVertexTest = 0; iVertexTest < 8; iVertexTest++) {
        if (afCubeValue[iVertexTest] <= fTargetValue) {
          iFlagIndex |= 1 << iVertexTest;
        }
      }
      return iFlagIndex;
    }

    public int getEdgeFlags(float fTargetValue) {
      //Find which vertices are inside of the surface and which are outside
      int iFlagIndex = 0;
      for (int iVertexTest = 0; iVertexTest < 8; iVertexTest++) {
        if (afCubeValue[iVertexTest] <= fTargetValue) {
          iFlagIndex |= 1 << iVertexTest;
        }
      }

      //Find which edges are intersected by the surface
      int iEdgeFlags = aiCubeEdgeFlags[iFlagIndex];

      //If the cube is entirely inside or outside of the surface, then there will be no intersections
      return iEdgeFlags;
    }
  }

  private class defaultFunction
      implements ImplicitFunctionInterface {
    @Override
    public double eval(float x, float y, float z) {
      return (float) Math.sqrt( (x - 1) * (x - 1) + y * y + z * z);
    }
  }
}
